## Replication material for the paper titled:
## Rural opposition to landscape change from solar farms: Explaining the diffusion of local ordinances against solar farms in South Korea
## Author: Inhwan Ko 
## Ph.D. Candidate in Department of Political Science, University of Washington
## Last Updated: Mar 23, 2023

### Load necessary packages
pkg_all <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
} # install & load at the same time

## List of packages
# for simcf and tile packages (made by Chris Adolph, University of Washington)
# access https://faculty.washington.edu/cadolph/index.php?page=60  

packages <- c("MASS","readxl","sf","sp","tmap",
              "raster","spdep","rgdal","rgeos",
              "ggplot2","ggpubr","tidyverse","survival", "spaMM",
              "simcf","tile","shiny","shinyjs","tmaptools") 

pkg_all(packages) # install & load

### Load data
load("localgovmap.RData")
load("data.RData")
load("adjacency_mat.RData")

extractCoords <- function(sp.df)
{
  results <- list()
  for(i in 1:length(sp.df@polygons[[1]]@Polygons))
  {
    results[[i]] <- sp.df@polygons[[1]]@Polygons[[i]]@coords
  }
  results <- Reduce(rbind, results)
  results
}

localgovmap_sp <- as(localgovmap, "Spatial")
centr <- gCentroid(localgovmap_sp, byid=TRUE)

### Piecewise event history analysis with spatial autocorrelation
data$panelperha <- (data$totalpanel / data$area) * 100
formula <- ~ panelperha + capaperpanel + logpopden + farmlandshare + netmig + margin + ghi
mdata <- extractdata(formula, data, extra=data, na.rm=T)

survdata <- Surv(mdata$year-1, mdata$year, mdata$reg)

# incorporate centroid coordinates into the mdata
coord <- as.data.frame(centr@coords)
coord$id <- c(1:89, 91:226)

for (i in c(1:89, 91:226)){
  mdata$x[mdata$id==i] <- coord$x[coord$id==i]
  mdata$y[mdata$id==i] <- coord$y[coord$id==i]
}

# Model 1: logit link with ICAR
formula_psf1 <- reg ~ panelperha*logpopden + capaperpanel*logpopden  + margin + farmlandshare + ghi + adjacency(1|id)
model_psf1 <- fitme(formula_psf1, mdata, family=binomial(link=logit), adjMatrix=E)
psf1_result <- summary(model_psf1)
psf1_result$beta_table %>% round(3)

# Model 2: logit link with Matern
formula_psf2 <- reg ~ panelperha*logpopden + capaperpanel*logpopden  + margin + farmlandshare + ghi + Matern(1|x+y)
model_psf2 <- fitme(formula_psf2, mdata, family=binomial(link=logit))
psf2_result <- summary(model_psf2)
psf2_result$beta_table %>% round(3)

# Model 3: cloglog link with ICAR
formula_psf3 <- reg ~ panelperha*logpopden + capaperpanel*logpopden  + margin + farmlandshare + ghi  + adjacency(1|id)
model_psf3 <- fitme(formula_psf3, mdata, family=binomial(link=cloglog), adjMatrix=E)
psf3_result <- summary(model_psf3)
psf3_result$beta_table %>% round(3)

# Model 4: cloglog link with Matern
formula_psf4 <- reg ~ panelperha*logpopden + capaperpanel*logpopden  + margin + farmlandshare + ghi  + Matern(1|x+y)
model_psf4 <- fitme(formula_psf4, mdata, family=binomial(link=cloglog))
psf4_result <- summary(model_psf4)
psf4_result$beta_table %>% round(3)

### model fit comparison: AIC
2*(2*10 - 2*model_psf1$APHLs$clik)
2*(2*11 - 2*model_psf2$APHLs$clik) 
2*(2*10 - 2*model_psf3$APHLs$clik) 
2*(2*11 - 2*model_psf4$APHLs$clik) # Model 4 best fit

summary(model_psf4)

source("theme_caviz.R") # custom-made theme by Brian Kai-Ping Leung, University of Washington
# tmaptools::palette_explorer() 


## Figure 1 ordinances contagion map, by year
rawdata <- read.csv("finaldata.csv")

# target year: input one year from 2015 to 2020
for (i in 2015:2020){
  target <- i
  subset <- rawdata %>% filter(year==target)
  figuredata <- left_join(localgovmap, subset, by="id")
  m <- figuredata %>% 
    tm_shape() +
    tm_fill("reg", palette=c("#F8F4F3", "#565555"), style = "cat", midpoint=NA,
            title="Local restriction", labels=c("Not adopted", "Adopted")) +
    tm_borders(alpha=.4)  +
    tm_layout(legend.position=c("right","bottom"))
  tmap_save(m, paste0("reg",target,".png"))
}


## Figure 2 interaction plot
# report 95% confidence interval
result <- summary(model_psf4)
covnames <- rownames(result$beta_table)
pe <- result$beta_table[,1]
vc <- vcov(model_psf2)
simbetas <- mvrnorm(10000,pe,vc)

coef <- colMeans(simbetas)
lwr <- matrixStats::colQuantiles(simbetas, probs=0.025)
upr <- matrixStats::colQuantiles(simbetas, probs=0.975)

hr <- exp(coef) %>% round(3)
hrlwr <- exp(lwr) %>% round(3)
hrupr <- exp(upr) %>% round(3)

data.frame(hr=hr, hrlwr=hrlwr, hrupr=hrupr)[-1,]

## set up counterfactuals
cfModel <- reg ~ panelperha*logpopden + capaperpanel*logpopden  + margin + farmlandshare + ghi
interactScen <- cfMake(cfModel, mdata, 5)

# scenario 1: higher solar farm density
interactScen <- cfChange(interactScen, scen=1,"panelperha", x=mean(mdata$panelperha)+1)
# scenario 2: higher solar farm density, higher population density
interactScen <- cfChange(interactScen, scen=2,"panelperha", x=mean(mdata$panelperha)+1)
interactScen <- cfChange(interactScen, scen=2,"logpopden", x=mean(mdata$logpopden)-1)
# scenario 3: lower solar farm density, lower population density
interactScen <- cfChange(interactScen, scen=3,"panelperha", x=mean(mdata$panelperha)+1)
interactScen <- cfChange(interactScen, scen=3,"logpopden", x=mean(mdata$logpopden)-2)
# scenario 2: higher solar farm density, much higher population density
interactScen <- cfChange(interactScen, scen=4,"panelperha", x=mean(mdata$panelperha)+1)
interactScen <- cfChange(interactScen, scen=4,"logpopden", x=mean(mdata$logpopden)+1)
# scenario 3: lower solar farm density, much lower population density
interactScen <- cfChange(interactScen, scen=5,"panelperha", x=mean(mdata$panelperha)+1)
interactScen <- cfChange(interactScen, scen=5,"logpopden", x=mean(mdata$logpopden)+2)

source("cloglogsimcf.R")
interactSims <- cloglogsimfd(interactScen, simbetas, ci=0.95)

labels <- c("+1 solar farm per ha",
            "+1 solar farm per ha \n & -1% population density",
            "+1 solar farm per ha \n & -2% population density",
            "+1 solar farm per ha \n & +1% population density",
            "+1 solar farm per ha \n & +2% population density")

col <- RColorBrewer::brewer.pal(5, "Set2")

interactTrace <- ropeladder(x=interactSims$pe,
                            lower=interactSims$lower,
                            upper=interactSims$upper,
                            labels=labels,
                            col=col[1:5],
                            size=0.8, lex=1.75, lineend="square", plot=1)

vertmarkFD <- linesTile(x=c(0,0), y=c(0,1), plot=1)
xat <- c(-0.05,0,0.05,0.10,0.15,0.20)
xlab <- c("-5%","0%","+5%","+10%","15%","20%")

file <- "figure2"
tile(interactTrace, vertmarkFD,
     xaxis=list(at=xat, labels=xlab),
     topaxis=list(add=TRUE, at=xat, labels=xlab),
     plottitle=list(labels="Figure 4. Differences in the likelihood of adopting the ordinance, \n Each scenario vs Baseline county"),
     xaxistitle=list(labels="Difference in the expected likelihood"),
     width=list(null=4),
     height=list(xaxistitle=3, plottitle=1),
     gridlines=list(type="xt"),
     output=list(file=file, width=8, height=8/1.618)   
)

## Figure 3 trend in deforestation
# data from http://www.e-platform.net/news/articleView.html?idxno=59289  (in Korean)
# Korea Forest Service reference to Member of Parliament Yoon Han-Hong

year <- 2006:2019
ha <- c(2,16,114,40,30,21,22,44,176,522,529,1435,2443,1024)

deforest <- data.frame(year=year, ha=ha)

deforest %>% ggplot(aes(year,ha, label=ha)) +
  geom_line(size=1.2, colour=col[3]) + geom_point() + 
  geom_text(hjust=0.6, vjust=-0.8, size=4) +
  xlab("Year") + ylab("Size of deforestation from solar farms (ha)") +
  ggtitle("Figure 2. Annual trend of deforestation from solar farms, 2006-2019") +
  theme(legend.position = "none", 
        plot.title=element_text(face="bold")) +
  theme_caviz_hgrid +
  theme(aspect.ratio = 1/1.1618, 
        axis.text.x = element_text(color = "black", size = 13),
        axis.ticks.x = element_line(size=1),
        axis.title.x.bottom = element_text(size=13),
        axis.title.y = element_text(size=13)) 

ggsave("figure3.pdf", width=10, height=10/1.518)

## Figure A1 Estimated correlation per distance between nodes
dd <- dist(mdata[,c("x","y")])
mm <- MaternCorr(dd, nu=16.666666667, rho=0.000126302) 

tibble(dd=dd/1000, mm=mm) %>% 
  ggplot() +
  geom_point(aes(dd,as.numeric(mm),col=col[1]),size=1) + 
  xlab("Distance between county centroids (km)") + ylab("Estimated correlation") +
  ggtitle("Figure 5. Estimated spatial autocorrelation against Euclidean distasnce between county centroids") + 
  theme_caviz_vgrid + 
  theme(legend.position = "none") +
  theme(aspect.ratio = 1/1.1618, 
        axis.text.x = element_text(color = "black", size = 10),
        axis.ticks.x = element_line(size=1))

ggsave("figureA1.pdf")

